% fig2a.m ——— 绘制 Fig.2(a1) 分岔图
clear; clc; close all;

%% —— 模型及仿真参数 —— 
a     = 0.5;       % 论文中 a'
b     = 0.215;     % 论文中 b'
beta  = 0.2;       % 论文中 β
A     = 0.95;      % 强迫幅值
omega = 0.61;      % 固定 ω = 0.61

T     = 200;       % 总仿真时间
dt    = 0.01;      % 积分步长
T0    = 100;       % 丢弃前 transient 时间

%% —— 参数扫描 —— 
alpha_list = linspace(1.6, 3.2, 150);   % α 从 1 到 4，共 200 点

% 用 cell 来存每个 α 对应的峰值数组
N_alpha = numel(alpha_list);
alpha_cell = cell(1, N_alpha);
ypeak_cell = cell(1, N_alpha);

X0 = [0.01; 0.02; 0.02];           % 初始条件
tspan = 0:dt:T;

for i = 1:N_alpha
    alpha = alpha_list(i);
    % 1. 用 RK4 积分系统
    X = rk4(@(t,X) rhs(t,X,a,b,beta,A,omega,alpha), tspan, X0);
    % 2. 丢弃 transient
    idx = tspan > T0;
    y   = X(2,idx);
    % 3. 提取峰值
    pks = findpeaks(y);            % pks 是列向量
    % 4. 存入 cell
    alpha_cell{i} = alpha * ones(1, numel(pks));
    ypeak_cell{i} = pks(:)';       % 转成行向量
end

% 合并所有点
alpha_all = [alpha_cell{:}];
ypeak_all = [ypeak_cell{:}];

%% —— 绘图 —— 
figure;
plot(alpha_all, ypeak_all, '.', 'MarkerSize',1);
xlabel('\alpha', 'FontSize',12);
ylabel('y_{peak}', 'FontSize',12);
title('分岔图：y_{peak} vs \alpha (\omega=0.61)', 'FontSize',14);
axis([min(alpha_list) max(alpha_list) min(ypeak_all) max(ypeak_all)]);
grid on;

